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We investigate the statistical mechanics of an inhomogeneous Coulomb fluid composed of charged 
particles with static polarizability. We derive the weak- and the strong-coupling approximations and 
evaluate the partition function in a planar dielectric slab geometry with charged boundaries. We 
investigate the density profiles and the disjoining pressure for both approximations. Comparison to 
the case of non-polarizable counterions shows that polarizability brings important differences in the 
counterion density distribution as well as the counterion mediated electrostatic interactions between 
charged dielectric interfaces. 



I. INTRODUCTION 

It has become clear by now that the standard Poisson-Boltzmann (PB) theory used to describe 
and understand the electrostatic interactions in colloidal systems has severe limitations and can 
sometimes give qualitatively unreliable if not outright wrong answers [l[ . There are several distinct 
reasons why the PB theory cannot describe some salient features of highly charged Coulomb systems. 

First and most notably, the PB theory is a mean-field theory, and thus completely misses the im- 
portant effects of ionic correlations that have recently been the focus of much research in Coulomb 
fluids 2]. The correlation effect, first observed in simulations Q, exposes the limitations of the 
mean-field ansatz in quite a drastic manner, since for highly charged systems the interactions me- 
diated by mobile ions between equally charged interfaces can become attractive. However, general 
theorems demand the interaction to be repulsive at the mean-field level 4-6]. Several lines of 
thought were spawned by simulations and converged into a paradigm shift that allowed for a simple 
conceptual understanding of why the mean-field picture breaks down for highly charged systems 
and how to formulate a theory that would be valid in these circumstances [7|. This paradigm 
shift led to a dichotomy between the weak and strong-coupling approaches that delimit the exact 
behavior of a Coulomb system at any value of electrostatic coupling Q . 

Another drawback of the PB theory is the physical model on which it is based - point charged 
particles - that neglects all ion-specific effects except for the ion valency. It is thus a one parameter 
theory where the ions differ only in the amount of charge they bear. One straightforward way 
to amend this drawback, sharing some of the conceptual simplicity with the original Poisson- 
Boltzmann formulation, is to take into account the excess static ionic polarizability of the ions 
proportional to the volume of the cavity created by the ion in the solvent. Static excess ionic 
polarizability is then a second parameter that differentiates between different but equally charged 
ionic species and thus presents an important step towards more civilized models of Coulomb fluids. 

Studies of the excess ionic polarizability have a venerable history and go all the way back to the 
classical book by Debye on polar molecules (see discussion on pages 111-115 in Ref. PJI]), where he 
already discussed cavities around ions having a different value of dielectric constant compared to 
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the surrounding solution. These cavities in fact represent excess polarization of the ions in aqueous 
solvent. Since due to saturation effects for most salts the interior dielectric constant should be taken 
much smaller than the aqueous one, the corresponding (static) dielectric constant of the salt should 
then be smaller than for pure solvent. This corresponds to negative excess ionic polarizability. 
While in Debye's analysis the effect scales as the volume of the ionic cavity, there are indications 
that for large enough solutes it should actually scale with the area of the cavity [13]. 

One of the moot points of Debye's analysis is exactly how to pick the right size of the cavity, an 
issue that has continued unabashed ever since [lij]. The changes in the effective dielectric constant 
of ionic solutions due to ionic polarizability was later picked up by Bikcrman [15] who, among other 
things, acknowledged that a realistic treatment of ions in aqueous solution should take their finite 
size into account (see also the discussion in [16]) as well as their excess polarizability. The effects 
of ionic polarizability and the associated dielectric decrement on the interactions between charged 
macromolecular surfaces in the presence of mobile counterions has been investigated in more recent 
times starting from the fundamental work of Netz [17| and continuing with a steady stream of works 

Some facets of the ionic and colloid polarizability were discussed starting from the weak-coupling 
level by generalising the zero Matsubara frequency van der Waals term and modifying the appro- 
priately formulated linearised Debye-Huckel theory [13, IH| ■ Levin and coworkers [Hi [2(| dealt with 
polarizability in the context of (ideally) polarizable ions in the vicinity of the dielectric interface. 
They also formulated a theory of monovalent and multivalent counterions in suspensions of polar- 
izable colloids or nanoparticles (2lj which in some respects complements our work where the mono 
or polyvalent counterions themselves are treated as polarizable. 

The main conceptual fulcrum of our present work is the dielectric decrement of ionic solutions 
that has been attributed to various sources, which underlie the changes in the dielectric response 
of the solution, but can be universally quantified by an excess ionic polarizability [H, It is 
proportional to the derivative of the (static) dielectric constant of a salt solution with respect to 
the concentration of the ions. Numerically this last coefficient, j3 9], turns out to be between 
-7M" 1 and -20 M" 1 for most of the common salts [2^ . 

Here we shall proceed with the analysis of effects of the excess static ionic polarizability of ions 
by formulating consistent weak- and strong-coupling approaches that will lead to a two parameter 
- charge and static excess polarizability - theory of a Coulomb fluid. We thus reformulate the 
basic model of a Coulomb fluid and investigate its consequences. This is accomplished by first 
incorporating the excess ionic polarizability effect in a consistent way into a microscopic model and 
then solving the corresponding theory at the mean-field weak-coupling level as well as at the strong 
coupling level. It further turns out that the radius of the ions (more precisely of their hydration shell 
or cavity) must be introduced, leading to a more civilized three parameters theory. The presented 
theory has thus a very broad parameter space that we can not analyze in complete detail. We point 
to some salient features and leave most of the details for future endeavor. 



II. MODEL 

We are interested in the behavior of mobile charges (counterions) immersed in a planar slab 
of thickness L filled by aqueous solvent of permittivity e w . The slab is assumed to be confined 
between two semi-infinite regions of permittivity e cxt that bear fixed charges of opposite sign to the 
sign of the mobile charges with surface charge density ctq- Counterions have a radius R, a charge 
e = qeo, where eo is the elementary charge of the electron and q is their valency, and an excess 
polarizability a. This latter quantity is defined precisely as the difference between the aqueous 
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FIG. 1: Polarizable counter-ions of excess polarizability a between two charged plates. The solvent in 
between has a permittivity e w , while the two semi-infinite regions > z > L have permittivity e cxt . The 
two surfaces at z — 0, L bear a surface charge of surface charge density a. R is the radius of the ions. 



solvent polarizability and the proper ionic polarizability, and may thus be negative as surmised by 
Debye [l2j. In fact experimentally this is the standard behavior observed for many salts, see Ref 
0] for details. We will denote the whole space as E and the volume of the slab as V . A schematic 
representation of the geometry of our model is given in Fig. [TJ 



A. Field-action 



The partition function for N countcrions is 




where f3 = (fc^T) -1 and d'x denotes the integration over the bounding surfaces dV. The standard 
field-theoretical representation of the Coulomb fluid partition function in terms of the fluctuating 
electrostatic potential has been used ;23j , properly extended by the fact that the counterion energy 
in an electrostatic field contains the point charge contribution ie<f>(xj) as well as the term due to 
its excess polarizability ^(V(j)(xj)) 2 . 

Note that in this general expression, the surface charge may not be uniform, although we will 
restrict ourselves to the case a(x) = a . 

The grand canonical partition function for a given fugacity A is then given by 

oc „ 

Z = xNz n = / exp(-/3S[0])[#], (2) 
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where the field- action S[<f>] is given by 

j dx — i(3 I a(x)(f>(x)d'x. 
' Jav 

This is the fundamental expression that we will evaluate; we will specifically concentrate on its 
dependence on the separation between charged plane-parallel boundaries. 



e(x)(V(f)(x)) 2 dx — A / exp 



(x) 



B. Dimensionless field-action 



The field-action can be rewritten in terms of dimensionless parameters. Of course this analysis 
holds only in 3D; in other dimensions, the characteristic lengths that define the dimensionless 



parameters would have to be defined differently [23[. The dimensionless form of the action itself 



suggests various approximations that allow an explicit and exact evaluation of the grand canonical 
partition function [7|. 

First, we recall the definition of the Bjerrum and Gouy-Chapman lengths, Ib = /3eo/4 7re weo and 
Iqc = l/27rqlB<Js, where <7o = e^us is chosen to be positive. The electrostatic "coupling constant" 
is then defined as the ratio [24| 



2irq i l B a s = q A S . (4) 



Above we specifically decomposed the coupling parameter into its q and as dependence. The 
dimensionless length, field, permittivity and surface charge can then be expressed as x = x/Iqq, 
<fi = (3ecj), e(x) — e(x)/e w , s(x) = —a(x)/a - One can also introduce a rescaled polarizability 
defined as 

a = 17> i — \2 a - ( 5 

{pqeolGcr 

Usually instead of using the excess polarizability one can use the dielectric decrement j3 in units of 
inverse Mole per liter [9[, defined as a = eo/3. Typically the dielectric decrement for various salts is 
negative. 

The dimensionless polarizability represents an additional independent parameter of the theory. 
Finally we define the dimensionless fugacity as A = 27tS^q C A. 

We can estimate the numerical values for all these parameters and obtain typical values for 
monovalent counterions that are of the order: Ib — lnm, S ~ 1, a ~ 10~ 2 and £ C xt — 5 x 10~ 2 . 

We can now derive the grand canonical partition function in the form 

Z = /exp(-^T) [#]) ( 6 ) 

where the field action can be obtained as 

S l<t>] = -tJ E < x )F<i>{x)fdx-^ j exp (-^(V^^)) 2 + i0O)) dx + ^-J^ s{x)cj>{x)d l x, (7) 

Here, in order not to proliferate the notation we simply renamed all the " ~ " quantities back to 
symbols, because in what follows we will work only with the dimensionless action. 
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This expression is then the point of departure for the evaluation of the free energy and pressure 
of the system. One should note here that the partition function ([6]) depends on two parameters: 
the coupling constant S as well as the dimensionless polarizability a, i.e. it is a two-parameter 
function. 



C. Density and electroneutrality 

In unscreened systems with long range Coulomb interactions the stability is insured only if the 
system as a whole is electroneutral. This is a particularity of long range interactions that becomes 
irrelevant for all finite range interaction potentials [25]. Special care then needs to be taken in order 
to stipulate this stability, that is given as a condition on the one-particle ionic density. The latter 
is defined by the operator 

n(x) = Acxp (-^{Vcp{x)f + i<p(xfj l v (x), (8) 

where ly(a?) is the indicator function of the volume V defined by J f(x)lv{x)dx — J y f(x)dx. 
We will also use the indicator function of the surface dV, given by J f(x)lgv(x)dx = J dv f{x)d'x. 

While the true density is actually given by g% , the above expression is easier to use in the mean 
field approximation since it does not involve S. We now impose average electroneutrality in the 
system by stipulating that 

/ (n(x))dx = / s(x)d'x. (9) 
J JdV 

This zero-moment gauge condition insures that the system remains stable for any configuration of 
the charges. Electroneutrality needs to be formulated as an additional condition on the density 
function only for unscreened interactions, sec 25] for details. 



D. Grand potential, free energy, and pressure 



The grand canonical thermodynamic potential is defined by 

J x = -lnZ x . (10) 

We write explicitly the dependence on A since it will feature prominently in our analysis. The 
fugacity is not a physical parameter, so the pressure should not depend on it. To solve this issue, 
we have to know how J\ depends on A; by differentiating ([6]), we get 

dX 2ttSA' [ ' 

where N has no subscript A because it does not depend on it as a consequence of electroneutrality 
Now, it is clear that the free energy defined by 

iV 

F x = J x + —\n\ (12) 
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does not depend on A, i.e. F\ = F, while In A is the chemical potential. This means that the free 
energy can be safely used to compute the pressure: 

dF 

Note that all the energies defined above are energies per unit area because of the transverse exten- 
sivity of our system. Furthermore the above pressure is in dimensionless units; the physical pressure 
is thus p = P/ipioc). 



III. WEAK COUPLING APPROXIMATION 



Depending on the strength of the Coulomb coupling as parameterized by the coupling constant H 
the grand canonical partition function exhibits two well defined limiting laws 0, l24j . For vanishing 
values of the coupling constant, 3 —¥ 0, the partition function can be well approximated by its 
saddle-point value and fluctuations around it. In fact the saddle-point is known to correspond 
exactly to the mean-field Poisson-Boltzmann expression while the Gaussian fluctuations around 
the mean-field correspond to the zero Matsubara frequency van der Waals or thermal Casimir 
interactions (26[. 

We will first derive the mean-field equations for our field-action, equivalent to those derived 
elsewhere [Tol ITlj . and then evaluate the Gaussian fluctuations around the mean-field and their 
dependence on the separation between the bounding surfaces. 



A. Mean-field 



We start with the general saddle-point equation satisfied at equilibrium 

~e(x) 



-V • 



+ an(x) 



V(f)(x) ) — i(n(x)) + is(x)lev(x) 



(14) 



The mean-field is more often written in terms of the (real) electrostatic potential proportional to 
ip = —i(f>, with the corresponding field-action S[ip], than in terms of the fluctuating potential <f>. 
For this new variable, the mean-field configuration is evaluated from the saddle-point condition 



SS[ip 



MF 



0. 



5xjj(x) 

The grand canonical potential is then approximated by 

S[4>mf] 



J ~ Jmf — 



(15) 



(16) 



From the saddle-point equation the mean-field equation can be rewritten in its Poisson-Boltzmann 
form as [H 



- + an } VipMF 



-n + slay, 



(17) 
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where the density is given by 

n = Aexp (^(VV>mf) 2 - <Amf) ly. (18) 

In these two equations, it is clear that the fugacity can be absorbed into the electrostatic potential: 
this change will modify the grand potential but not the free energy. We can thus assume A = 1 for 
the mean-field as well as for fluctuations around it. 



B. Pressure in the plane-parallel geometry 



In ID, which is also the case of two charged plane parallel surfaces since the mean potential 
depends only on the transverse coordinate z, the Poisson-Boltzmann equation has the form 



with 



(| + an(z)J iPmf{z)' = -n(z) + sl dv (z) 
i(z) = Aexp (^iPmf(z)' 2 - ipMF(z)j lv(z). 



(19) 



(20) 



We used the notation f'(z) = ir{z)- In this case it can be shown that the pressure in the system 
is a constant given by the contact value theorem Q 



P = 



1 



2ttS 



n (z) - (| + an(z)j ijy MF (z? 



= const. 



(21) 



It can be easily checked that this quantity is actually equal to the pressure obtained equivalently 
by the standard thermodynamic definition P = — = dShW] ^ a ^ ove f orm f or the interaction 

pressure contains an osmotic van't Hoff term, the first one in Eq. [51] that contains the effects of the 
polarizability implicitly, i.e. through the variation of the density profile on the polarizability, and 
a Maxwell stress term, the second one in Eq. [5TJ that contains the polarizability effects explicitly. 



C. Second order fluctuations correction 



The grand potential can be computed to the next order by taking into account fluctuations 
around the mean-field solution. This is done by expanding S around 0mf = *"0mf to the second 
order, obtaining 



S[(pMF + 0}= <S[V>mf] + 



S 2 S 



S(f>(x)S(f)(y) 



[4>MF}0(x)0(y)dxdy = 5[^mf] + S^[0}. 



In this case the grand potential is given by 

T (i) _ S[i/j M f] _ , _ (2 ) _ S[iPmf] 



J ~ J, 



MF 



In 



exp 



SW[0] 



[d0] 



(22) 



(23) 
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where is the contribution of the fluctuations to the partition function. The effective action for 
the fluctuations [6] is straightforward to compute, yielding 



SW[6] = i- J (| + an) (V0) 2 + na 2 (V^ M F • V0) s 



-V • (eV^mf) 



slay 9 



, (24) 



where we have used the mean- field Eq. (|T7|) . 

Note that eVt/'mf is not continuous, and thus leads to a surface term. As expected, this ac- 
tion does not depend on the fugacity but on the mean-field density, which is the only physically 
meaningful quantity. 

In our model we consider parallel plates of constant surface charge, the mean field problem is 
then one dimensional. We can therefore split the coordinates in a one dimensional coordinate z 
perpendicular to the plates, and a two dimensional one parallel to the plates: x = (z,r). 



D. Pressure 



Fourier transforming the fluctuations in the direction parallel to the plates we obtain 

dk 



9(z,r) = / exp(ifc • r)6(z, k) 



(2tt) S 



(25) 



where 9{z, —k) = 6{z, k)* because 9 is real. This decomposition furthermore allows us to write the 
fluctuations action (pM|) as 



= s 



:(2) 



9(;k) 



dk 



(26) 



where the one dimensional action is 



S k 2) [9] 



1 

47T 



+ an + a nip' 



MF 



— [9{Qf + e{Lf 



_ c(2) , o(2) 
— D fc,b "+" °fc,s- 



(#mf)' + [- + an) k 



(27) 



(2) (2) 

This action thus has a bulk part Sj, L and a surface part S k ' . The surface action actually contains 
another term due to the fact that stPmf ^ s discontinuous across the bounding surfaces, so that finally 



S k 2 l[9} = ^(9(0) 2 + 9(L) 2 ), 



c 



where 



2tt' 



(28) 



with the notation [/(x)]^ = f{x%) — f{x\). We also used the symmetry zni-zof our system. 
The partition function for the fluctuations can be written as a product of path-integrals, 



z(2) = n/ ex p 

k 



Si 2) [6] 



(29) 



9 



These path- integrals are computed in appendix [XJ leading to 



z k = ex P Hr 



2 / 



, (30) 

C+e^tk/iTT 1 ^ 



\ a f fc (0,£)+ c+£ °g fc/4 " -6 fc (0,L) 



where the functions 6 fc and af are defined in the appendix. 

The total free energy of the mean field configuration and fluctuations around it is then obtained 

as 

F& = ^-JLf o °° m (3i) 

We note here that the structure of the free energy Fj\l does not look like a mean-field term 
independent of the counterion polarizability plus a zero frequency van der Waals term that stems 
from the polarizability of the counterions. Though this kind of decomposition is sometimes assumed 
in the literature [27L [28[ , it clearly does not correspond to the weak-coupling approximation. 

We can see numerically that the integral over the transverse Fourier modes in Eq. (|3"Tj) diverges; 
this comes from our model of point-like dipoles. Taking into account the size R of the polarizable 
ions (more precisely, R is the charge of their hydration shell), the integral is regularized by the 
dimensionless cut-off 

^max — 7> " (^) 
K 

Physically the cut-off arises because electric fields which fluctuate on length scales shorter than the 
polarizable ion cannot polarize it. The interaction pressure on this level of approximation is then 
obtained by taking into account Eq. leading to 

p(i) = _^M£_ (33) 

The results for the fluctuations-corrected interaction pressure from Eq. [33] on the weak coupling 
approximation level are shown on Fig. [5] for S = 1, e cxt = 0.05 and R = 1, for various values of 
the counterion polarizability a. The fluctuations correction in P*- 1 ) is quite small compared to the 
mean-field value, but can become substantial as the polarizability a decreases, i.e. becomes more 
negative. This correction reduces the interaction pressure between the surfaces. This indicates 
that ions with nominally equal charge (of equal valency) but differing in the polarizability will 
mediate markedly different interactions when confined between charged dielectric interfaces even 
at the weak-coupling level. 



E. Density 

We now consider the ion density by taking into account the mean field solution as well as the 
fluctuations around the mean field. From (1181) . the ion density is 

Pl (x) = (exp(-|(V^)) 2 + z0( : r))) i (34) 
= (27ra)" 3/2 J exp (^-^) ( ex P (iP ' V0(cc) + i(f>(x))) 1 dp, (35) 
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FIG. 2: Pressure P m {L) fr om Eq. [33] as a function of the plate separation L for e C xt = 0.05, S = 1, and 
R = 1. The dashed lines are the mean field result, the solid lines include the fluctuations. Left: for a = —0.1 
and large plate separation, the difference is barely distinguishable. Right: for small plate separation L and 
various values of the polarizability a. The effect of fluctuations can be quite important for large counterion 
polarizabilities. 



where we used a Hubbard-Stratonovitch transformation to obtain the last expression. In this way we 
have only terms linear in <fi in the exponential. The subscript 1 denotes that we take into account the 
first order of the fluctuations. We notice that the mean- field equation for electroneutrality should 
also hold on average at equilibrium, so that p\ will satisfy electroneutrality. As a consequence, we 
only need p\ up to a multiplicative constant, and this constant will be set by electroneutrality. 

The interpretation of the above formula is that the local ion density is the average over a fluctu- 
ating dipolar moment vector of a Coulomb fluid characterized by ions with a charge and a dipolar 
moment. We then decompose <f> into a mean- field term plus Gaussian fluctuations 

4> = #MF + (36) 

obtaining 

pi(x) = (27ra)" 3/2 J exp - p ■ Vip MF {x) - ^uf(x) \ (exp (ip ■ V0(sb) + i9{x))) 1 dp. (37) 

The average is now easy to compute, 

(exp (ip ■ V8(x) + i0(x))) 1 = exp f-i <(p ■ V0(x) + 9(x)) 2 ))j , (38) 

and then, introducing the correlator of the fluctuations 

G(x,x') = (6(x)6(x')) u (39) 

we can write it as 

(exp(ip- V9(x) + i6(x))) 1 = exp (-^p T \7\7' T G{x, x)p - ^G(x, x) - -p-VG(x,x)\ . (40) 
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We used the notation V for the gradient with respect to the first variable of G(x,x'), V for the 
second variable, and V for the sum of the two gradients. We can now insert this expression into 
(137|) , and remain with a Gaussian integral 



pi(x) = (2ira) 3/2 J exp (^-^p T a 1 l(x)p - p ■ Vipi(x) - ipi(x)j dp, 



(41) 



where we introduced a renormalized polarizability (which is now a position-dependent matrix) and 
a renormalized field 

a^(x) = a' 1 + VV' T G{x,x), (42) 
^(x) = ip M F(x) + ~G(x,x). (43) 

Performing the integral gives 



Pl+ (x) = ^dct (^^) cxp QlVTM^fai+NVVM*) - M*)) ■ (44) 

The index "+" means that our computation works only for a > 0. In the more common case where 
a < 0, the computation is the same up to some factors of i, and we get 

olJ1(x) = |a| _1 - W T G{x,x), (45) 

and 



Pi-{x) = ^dct (^^j cxp ^{VM*)} T <Xi-(x)VMx) ~ M*)) • (46) 

Now we have to compute G(x,x) and W' T G(x,x) at each point. To do this, we will use the 
same technique we used to compute the pressure: we Fourier transform the fluctuations in the 
direction parallel to the plates and use the Pauli-van Vleck formula. 

Since the fluctuations action (|24|) is a sum over different transversal modes, two modes with 
different wave vectors are uncorrelated and we can write the correlator as an integral over the 
modes. 



r dh 

G(x,x') = J e^(ik-[r-r'])G k {z,z')—^, (47) 



(27 

where Gk(z,z') is the one dimensional correlator for the action in Eq. (l27l) . More precisely, we 
need Gk{z,z) as well as dd'Gk(z, z). These functions are computed in appendix |B| Then we can 
write 



/dk 1 f°° 

° k ^ ^ (27r]2 = 2tt J G k {z,z)kdk, 



(48) 



and the matrix 



VV- G <*,*> - j a( „„^ I - ^ jf (T I) ,49, 
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FIG. 3: Counterions density profile within the slab - dependence on the coupling constant S. Left: Weak 
coupling density close to the left electrode taking into account the fluctuations around the mean field as 
a function of the position within the slab (z 6 [0, ■§■]), for R — 0.3, a — —0.3, e cx t = 0.05 and S = 0.3 
(dashed line) or S = 1 (solid line). The mean-field density itself is presented by the dotted line. Right: 
Strong coupling density as a function of the position for R — 2, a — —0.01, e cx t = 0.05 and S = 10 (solid 
line) or 5 = 50 (dashed line). 



where I2 is the two dimensional identity matrix. 

In conclusion, we have the algorithm to compute the pressure and the density: for each mode, 
we integrate the Pauli-van Vleck formula and then compute its contributions to G(x, x) and 
VV /T G(a;, x) and add them to the contributions of the previous modes. Finally, we compute 
the new density at each point and renormalise it using electroneutrality. 

The mean field and first order densities can be compared on Fig. [3] (left). First of all, we observe 
that the effect of the fluctuations is small and depends on the values of the parameters. For higher 
2, e.g. 5 = 1 on the figure, the ions get preferentially included in the region close to the dielectric 
boundaries. This is not a mean-field effect since the mean-field density does not depend on the 
coupling parameter. The a dependence of the counterion density in the slab is shown on Fig. 2] 
(left). The mean field density depends strongly on a [2j, |9(, and this dependence remains after one 
adds the fluctuation contribution. The inset shows that the deviation from the mean-field density 
increases with a: the effect of the fluctuations is enhanced by the polarizability. 

As a conclusion, the counterions are attracted by the boundaries, and most of this effect has a 
mean-field nature. The polarizability tends to lower this attraction at the mean-field level, but to 
increase it at the fluctuations level. 



IV. STRONG COUPLING 

The strong coupling approximation 0, [24| is obtained in the limit of asymptotically large coupling 
parameter, E — > 00. In this limit it turns out that the statistical mechanical description of the 
system is equivalent to a properly normalized one-body description. This means that we can treat 
the system as composed of bounding surfaces and a single polarizable charge between them. We 
will first derive the strong coupling form for the partition function, equivalent to the first order 
virial expansion, and then evaluate the density profile and the interaction pressure. 
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FIG. 4: Counterfoils density profile within the slab - dependence on the polarizability a. Left: Weak 
coupling density close to the left electrode taking into account the fluctuations around the mean field as a 
function of the position within the slab (z 6 [0, ■§■]), for R = 0.3, £ cx t = 0.05, H = 0.3 and a — —0.3 (solid 
line), a = —0.1 (dashed line) or a — — 10~ 6 (dotted line). Inset: Deviation from the mean-field density. 
Right: Strong coupling density as a function of the position within the dielectric slab for R = 1, e cxt = 0.05, 
S = 10 and a = -0.05 (solid line), a = -0.01 (dashed line) or a = -10" 6 (dotted line). 



A. Formulation 

The strong coupling limit formally corresponds to the A < 1 limit. To the lowest non-trivial 
order in A, the partition function is given by 



X 



2vrS 



-Z x =Z [1 



X 

2^S' 



(50) 



We are thus interested in the evaluation of 



Z = [d<t>] exp 



2ttS 



- / £(x)(V(t>{x)fdx + i / s(x)4>{x)d' x 
4 J Jdv 



and Z\ — J zi(xo)dxQ, with 



(51) 



zi{x ) 



/[#] 



exp 



2vrS 



s{x){\Jcj)(x)) 2 dx + i / s{x)4>{x)d' x 

4 J JdV 



-(V0(^ o )) 2 -i<p(x ) 



(52) 



The quantity Xz\(xq) / Z ~ \z\(xq) / 'Zq is the ionic density at Xq- The total number of ions thus 
follows by stipulating that XU = N, so that A can be tuned to satisfy electroneutrality. 

As in the mean-field approximation, we will be interested in the density and the pressure. For 
the density, we will be specifically interested in the Xq dependent part of zi(xq), whereas for the 
pressure we need the L dependent part of Zq and Z%. In this sense the density is easier to compute, 
so that we address this question first. 
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B. Density 

We introduce an auxiliary vector p together with a Hubbard-Stratonovich decomposition and 
perform the integration over cj> to write down Eq. [52] as 



Zl (x ) = ( 27 rar 3 / 2 det(-g^ 



2\ -1/2 



x J dpexp^-^-i^-^ J s(x)cp(x)d'x+p-V<f>(x ) + cj>(x )^ ^ ,(53) 

where (• ■ -)o denotes the Gaussian average over with the action So[<p] = / s(x)(V<fi(x)) 2 dx. 
In this way the average can be written in terms of the correlator, G(x, x') = (<f>(x)<f>(x'))o, as 

/ K^W^+P'WW + ^wj | =p T A(x )p + 2p-B{x ) + C(x ), (54) 

where 

A(x ) = VV' T G(x ,x ), (55) 

(56) 



B(x Q ) = V (g{x ,x ) ~ ^A=r J s(x)G(x ,x)d'x) , 



C(x ) 



G(xo,Xq) [ s(x)G(xo,x)d'x + — / s(x)s(x')G(x,x')d'xd'x' 

C"(a:o) + ,„ * [ s(x)s(x')G(x,x')d'xd'x', 
(2n^) z J 



(57) 



and V and V denote respectively the gradient with respect to the first and second variable. We 
can now perform the explicit integration over p, obtaining 



z 1 (x ) = Z Q dct(l + aA(x )) ^xexp^-BixoY ^- + A(x )j B(x ) ^ I . (58) 

where 

= dct (--^) x ex P (~ 2 (2t!") 2 / s ( x ) s ( x ') G ( x i x') d ' a;d ' aj ') ■ ( 59 ) 

As we noted in the weak-coupling treatment, the Hubbard-Stratonovitch transform depends on 
the sign of a. However, it is easy to see here that the final expressions (|55ti58|) remain the same if 
a is negative. 

We should mention that a problem arises in Eq. (|58|l if an eigenvalue of 1 + aA(xg) is negative. 
This happens notably if a is too negative, leading to a negative effective permittivity of the hydration 
shell of the ion, and thus to an instability for the field. The problematic value of a therefore strongly 
depends on the radius of the hydration shell. 

In order to be more explicit, we need the expression for the correlator. Again, we Fourier 
transform the field in the direction parallel to the plates: 

dk 

(2V) 



(z,r)= I cxp(zfc ■ r)4>(z,k) ^, d _ 1 , (60) 
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and the correlator for the k mode is relatively easy to determine and is given in 29] . To make the 
symmetry z — > L — z more explicit, we switch to coordinates where the plates are located at — L/2 
and L/2; in this case the correlator is given by 



G k (z,z') =4ttS 



exp(— k\z — z'\) cosh(fc(,z + z')) + A exp(— kL) cosh(fc(z — z')) 



2k 



where 



A = 



A- 1 exp(fcL) - A exp(-fcL) 



1 — £oxt 



1 + £oxt ' 



(61) 



(62) 



We will write A(xo), B(xq) and C'{xq) using this expression for the correlator. Divergences may 
appear, but for the density itself we can drop (almost) all the :co-independent terms. In fact we 
find 



A( Xo ) = ± 



dd' 




G k {z ,z )kdk, 



(63) 



where we need a cut-off as in the weak-coupling limit, and 

(k) 



(fc) — k arctan 



dd'G k {z ,z ) = 4S 
where q max (k) is defined by (|B9[) . q 
B(x ) 



c(fc) 5 



k 
k 2 



,_=/.2 cosh(2fcz )- Aexp(-fcL) 

" A- 1 exp(fcL)- Aexp(-fcL)' K ' 

Then, for B(xq), we will keep 



k 2 sinh(2fcz ) 



A" 1 exp(fcL) - Aexp(-fcL) 



dk. 



and we can drop the second term in C"(xq), 



C'(x ) = 2~ 



cosh(2fczo) 



A" 1 exp(fci) - A exp(-fci) 



k dk. 



(65) 



(66) 



The result is shown in Fig. [3] (right), where we used electroneutrality to normalize the strong 
coupling result, defined as we have seen up to a constant. We see that the counterions are completely 
excluded from the region close to the interfaces and pushed towards the middle of the dielectric 
slab. This effect increases with the coupling parameter S. It appears on Fig. @]that, in opposition 
to the weak coupling limit, the dependence on the polarizability is weak. Fig. [5] shows that the 
strong coupling density is ruled by the images: without them, the density would be constant within 



the slab 



24, 3C 



As a conclusion, the polarizability has a small effect at the strong coupling level. 



C. Pressure 



Using its definition p2|) and the condition XU = N, we can write the L dependent part of the 
free energy, up to the first order in A, as 

N N N 

F = J + ^lnA = -InZ -lnZ7~ -1nZ -InU. (67) 

2ttzl 2tt^ 2ttzl 



1G 



Of course we need to take into account the zero-moment gauge condition (electroneutrality, EqJSJ) 
when evaluating the above expression, which to the lowest order eliminates the S dependence. 
Let us first compute the L dependent part of Jq = — In Z , using (|59p we get 



Jo =2 In 



det 



47rS 



2(2tt 



^— T [ s(x)s(x')G(x 7 x')d'xd'x'. 



(68) 



The first term is the thermal Casimir fluctuation free energy, and the second the electrostatic 
interaction between the plates. Using the decomposition of the correlator in orthogonal modes, we 
can write 



s(x)s(x')G(x, x')d'xd'x' = 2 (G o (0, 0) + G o (0, L)) 



(69) 



where we have taken J dr = 1 to have the energy per unit area. Since the orthogonal mode k = 
is ill-defined in (|6"T1) . we can derive with respect to L before taking the limit k — > 0. We get 
^-G fe (0,0) -> and -£-G k (0,L) -> -2nE, so that we can replace G o (0, L) = -2vrSL. We thus 

get the final expression for the grand potential, 

Jo = — / In (1 - A 2 exp(-2kL)) kdk - 7 —. (70) 

47T J 27T^ 

Within the strong-coupling virial expansion this term corresponds to the free energy of the system 
without any ion. 

Now we need to compute U, remembering that we only need the L dependent terms in In U. 
Starting from expression (|58[) . we get 

zi(x ) = Z exp(-W(x )), (71) 

where W(xq) is an effective one-body potential given by 



W(x ) = ~B(x ) 



A(x ) B(x ) 



C'(x ) 



iTrlog(l + aA(a;o)). 



(72) 



The Zq term and the last term in the exponent of the above equation can be rearranged and 
interpreted in the following way: keeping only terms proportional to (V</>) 2 in the exponential of 
(l52"l), we can show that 



det 



eV 2 
47rS 



-1/2 



det (1 + aA(x Q )) 



-1/2 _ 



= det -V 



£ ( X ) cv \ 

— + ao(x — xq) 



AttE 



-1/2 



(73) 



This means that these two functional determinants represent the thermal Casimir partition function 
for a system composed of a finite extension dielectric slab, two semi-infinite dielectric regions outside 
of it and a single polarizable ion within the slab. The delta function in the expression for the effective 
dielectric response function on the r.h.s. of (1731) needs to be regularized to avoid a divergence in the 
case of a point ion. It is clear that the last term in the effective one-body potential (|T2")) describes 
the thermal Casimir or zero-frequency van der Waals interaction between the polarizable particle 
and the dielectric interfaces in the system. It is given explicitly by 



iTrlog(l + aA(x )) = ^Trlog (l + aVV' T G(;co, sc„)) ~ ±aTr [VV' T G(cc , ^o 



(74) 
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In the asymptotic regime of large Xq we obtain the scaling x^ 3 which corresponds to the zero- 
frequency van der Waals interaction between the polarizable particle and a single dielectric discon- 
tinuity [3l|. Our results are thus completely consistent with everything else we know about the 
polarizable particles and their zero-frequency van der Waals interactions with dielectric discontinu- 
ities d3. 

At the end of the calculation we obtain for the L-dependent interaction free energy (|67|) first the 
usual extensive term, A, giving rise to an attractive force between the plates, independent on L 
while the other terms can not be evaluated analytically but are easily calculated numerically: the 
computation of A follows from (|63|) . B is obtained from (|65|) and finally C'(xq) is obtained in the 
form 

where we differentiated with respect to L before taking the limit k — > in order to get the last 
term that represents the attraction between the ion and each of the bounding dielectric surfaces of 
the slab. Finally U is given by integrating (f72|) . 

The pressure obtained from the interaction free energy is shown on Fig. [5] for different values of 
a and R. The dependance on the polarizability is very weak and non-monotonous, contrarily to 
what we observed in the weak coupling limit. The results are however very sensitive to the size of 
the ions, especially for large ions. 



V. DISCUSSION AND CONCLUSIONS 



In this paper we have formulated a theory of Coulomb fluids that, apart from the charge of 
the mobile counterions, includes also their static excess polarizability. This leads to a possibility 
of ion specific effects even for ions with nominally equal valency |9|. Instead of starting from 
the phenomenological description of the ionic effects on the local dielectric function, an endeavor 
pursued in Ref. [g.[Io| , we rather implemented the effect of ionic polarizability at the level of the field 
action deriving the appropriate field-theoretic representation of the model. Though this variation 
in the approach results in the same form of the model at the mean-field level, the formulation 
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FIG. 6: Effect of the outer permittivity on the counterion density distribution. Left: Weak coupling density 
with a = —0.3, R = 0.3, 3 = 0.3 and e cxt = 0.05 (solid line) or e cx t = 1 (dashed line). Right: Strong 
coupling density with a — —0.01, R — 2, S = 10 and e ex t = 0.05 (solid line) or e cx t = 1 (dashed line). 





FIG. 7: Effect of the outer permittivity on the interaction pressure. Left: Weak coupling pressure with 
a — —0.1, R = 1, S = 0.3 and e oxt = 0.05 (solid line) or e cxt = 1 (dashed line). Right: Strong coupling 
pressure with a = —0.01, R = 2, S = 10 and £ cxt = 0.05 (solid line) or e cx t = 1 (dashed line). 



presented is eventually more general and suitable for further analysis and implementation of the 
weak and the strong-coupling asymptotic limits. 

After formulating the model and casting it into a field-theoretic form, we derive the pressure and 
the ionic density in the mean-field level approximation - corresponding to the saddle-point of the 
field-theoretic action. We then add the effects of Gaussian fluctuations of the local electrostatic 
potential around the mean field. This constitutes the weak-coupling approximation of the complete 
field theory. 

The effect of the fluctuations around the mean-field saddle point is found to be rather small. In 
the pressure itself it is barely discernible, see Fig. [3J but it does become stronger as the polarizability 
of the ions is increased, while the density profile shows an effect only very close to the boundaries of 
the system where the ionic density is enhanced, see Fig. [3J depending on the coupling parameter. 
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This modification of the ionic density in the region close to the dielectric boundaries of the system 
is partly due to the image effects [29l . |30| and partly due to the ionic polarizability. 

We then formulated a full strong-coupling theory which formally corresponds to a single-particle 
level description and derived its consequences in detail. In the strong-coupling limit the ions are 
expelled from the vicinity of the dielectric boundaries, see Fig. [3l The origin of this effect lies in 
the dielectric i mag e interactions that lead to a vicinal exclusion of the ions close to the dielectric 
discontinuities [29L l3Cj. 

The results derived here, both for the mean-field plus fluctuations and for the strong coupling 
regime, exhibit a dependence on the ion polarizability as well as on the size of the ions, Fig. [2] 
and The effect of ionic polarizability on the interaction pressure is connected partly with the 
changes in the density profile leading to the changes in the osmotic van't Hoff component of the 
interactions pressure, and partly with their contribution to the Maxwell stress term in Eq. 1211 The 
ionic size dependence comes from divergences naturally present for point dipoles. We have to stress 
here that what we refer to as the size of the ions is actually the size of the ionic cavity in the solvent 
which includes also their hydration shell Q • Despite the fact that the field theory arising from our 
model is a priori independent of the ionic size, we see that if one leaves the domain of mean field 
theory, by either taking into account fluctuations or going to the strong coupling limit, calculated 
thermodynamic quantities exhibit ultra-violet, or short distance, divergences. These divergences 
are associated with the inclusion of ionic polarizability as they do not arise in strong coupling or in 
the mean field fluctuations for non-polarizable ions. We have argued that the length scale used to 
cut-off ultra-violet divergences is thus the size of the polarizable molecules. Our results are thus in 
line with Bikerman 15] who long ago argued for the role of the ionic size. The effects of the ion size 
on the interaction pressure are shown in Fig. [5] The size dependence mediated by the polarizability 
of the ions has nothing to do with steric effects and has not been seen before for non polarizable 
ions or for polarizable ions on the mean-field level. 

The effect of dielectric images, i.e. of the outer permittivity, is shown in the case of density on 
Fig. |5] and in the case of pressure on Fig. [7] The weak coupling limit is only weakly affected by 
the images, which is to be expected since this regime is dominated by the mean field that does not 
depend on e oxt [2{|. On the other hand, the strong coupling limit is strongly affected, this time 
because images add a non negligible term to the correlator (I5T1) . 

Due to the polarizability of the ions, it is also clear that our two approximations break down if 
the parameters are too extreme, but for different reasons. This is easy to analyze for the strong 
coupling result (|58[) : here extreme parameter values correspond either to ions which are too small 
or have too high a (negative) polarizability. In this case, the effective permittivity around the ion 
may turn negative, leading to a field instability that shows up in the partition function. For the 
weak coupling limit, on the other hand, if the dielectric function 1 + an(z) becomes negative on 
the mean field level, a divergence appears for the fluctuations about the mean field and the system 
becomes unstable. This can be interpreted in the sense that the validity of the strong coupling vs. 
weak coupling description no longer depends on a single coupling parameter, but actually on three 
parameters. More work would thus be needed to explore different regions of the parameters space 
and assess the validity of the WC-SC dichotomy in each of them. Our present work can only be 
seen as a first step towards this complicated endeavor. 

One general conclusion stemming from the present work is that the contribution of polarizable 
counterions to the total partition function is in general non-additive, contrary to what is sometimes 
assumed [13, HH [32[. It is in fact highly non-additive at the weak coupling level, whereas it can 
sometimes be reduced to an additive contribution in the free energy at the strong coupling level, only 
if the polarizability is large enough. Simply adding a van der Waals ion-polarizability dependent 
contribution to the electrostatic potential of mean force is simply wrong. 
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A final note is in order about the possible computational verifications of our analytical calculations 
via coarse grained simulations that we did not attempt to see through in this work. As polarizability 
belongs to non-pairwise additive effects the simulation of the present model presents a considerable 
challenge. One would need to include the image interactions as well as the polarizability couplings 
to all orders which would appear to be no small accomplishment. Until such time when these type of 
simulations are actually performed our analytical calculations will remain the sole means to assess 
the consequences of our model of Coulomb fluids. 
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Appendix A: Fluctuations partition function 



In this appendix we will compute the path-integral appearing in (I29|) , 



2f= /expf-^W (AD 



We introduce the propagator 



f 9{*>)=Bx ( sfi\e\\ 

K(9 ,d 1 ;z,z';k) = exp *±LL\[d0), (A2) 



8(z)=8 

where, implicitly, the action is taken only over [z, z'\. 

We denote \w_K ext {9o, @i] I] k) the propagator for a mode k on a length I in the external medium; 
it is given by 33 



- V 8^SmM 6XP l- S^taW ^ + + ^Mkf) 6001 )- (A3) 

We also denote K^(9 0l 9i; z, z'; k) the propagator between z and z' in the ionic solution with the 
field tp (and the ionic density n contained implicitly in the field). The path- integral (jAll) is thus 

Z< 2) = Km J K cxt {9 ,9 1 -J;k)K^ MF {9 1 ,9 2 ;0,L;k)K cxt (9 2 , 9^1^)0^^-^ [9\ + 9%]^J\dti j , 

C9 2 \ „ ,„ , _ / C6» 2N 



= \hn(l\K ext (l;k)exp^- — J ^ MP (0, L; fe) exp ) K cxt (l;k)\l), (A4) 

where we introduce a matrix notation for the propagator. We have to integrate over the possible 
outer values of the field, so we need 

K ext (9 o ,9 i; l;k)d0 o = . L== exp ( - g «t* ta ^(*0 (A . (A5 ) 

(2) 

We can see immediately that the limit I — > oo is not well defined here: we find Z^> = 0. bmce 

(2) 

we need only the L dependence of Z k up to a multiplicative term, we can remove the term 
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l/ycosh(fc/) in the equation above. Thus the path- integral becomes 



Z< 2) = lim 

I— »oo 



cosh(fcZ) ^1 
= (ext, fc|iC^ MF (0, L; fe)|ext, k), 
where we have used 



CO 2 \ ( CO 2 

K ext (l;k)exTp ( — — ) i ; T l/ , MF (0,i;fc)exp ( — — ) K ext (l;k) 



2S 



2S 



| ext, fe) = exp 



C + e cxt k/Aw 
2S 



To evaluate the propagator between the plates, we write the bulk action 



with 



Sj?j[8] _ 1 
3 2 

.4 



[A(z)0'(z) 2 + B k (z)9(zf] dz, 



1 M 2 ,« 

— I — h an + a 



2vrS \2 



MF I ' 



2vrS 



MF 



Then the propagator is of the form [35 



K(e ,e 1 ;z,z';k) = J bk{ *'_ Z ' ) e>. V 



2tt 



— + an ) fc^ 



- "i^lel + b k ( z , z') 



It is easy to show that a k , a k , and 6 obey the following composition rules: 



a k (z,z' + () — a k (z,z') 



b k (z,z') 2 



a k {z,z') + a k {z',z' + QY 



a k (z,z' + () = a k (z',z' + C) 
b k (z,z' + 



b k (z',z' + C) 2 



a k (z,z')+a k (z',z' + ()'' 
b k (z,z')b k (z',z' + C) 



af(z,z') + af(z',z' + C) 



1 



(A7) 



(A8) 

(A9) 

(A10) 
(All) 

(A12) 

(A13) 
(A14) 
(A15) 



On the other hand, one can show that on a small interval [z, z + Q, where A and B k are almost 
constant, they are given by 



i k (z,Z + Q - 



^A{z)B k {z) 
tanh(w fc (z)C) 



a k (z,z + C) = 
1 ' s ' tanh(w fc (2;)C) 

b k (z,z + = 



sinh(cjfc(z)C) ' 



(A16) 
(A17) 
(A18) 
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where u>k(z) = y/Bk(z)/A(z). From eqs. (|A13IIA18|) we can show that a k , af and b h satisfy 



db k , _ _a k {z,z')b k {z,z') 
t (z,z ) — 



dz 

oz 
da k 



A(z') 
b k {z,z') 2 

W) 1 



= B k (z')- 



a k (z, z') 2 
A(z') 



with the initial condition 



b k (z,z') ~ a k (z,z') <~ a^(z,z') > 



A(z) 



z'^>z Z — Z 



(A19) 
(A20) 
(A21) 

(A22) 



Eqs. (|A12|) and (|A19flA21|) are the Pauli-van Vleck formula. We will however use eqs. (|A13tiA18|l 
for the numerical integration. 
Using (|A12|) in (|A"6|) . we get 



r(2) 



27rfo fc (0,L) 



\ k(o,i)+ 



C + £cxtfc/47T 



6 fe (0,L) 2 



(A23) 



However, this expression leads to a pressure that contains a constant term (i.e. independent on 
L), that comes from the fact that in our computation the total volume of the space depends on L. 
This term is thus not physical, and should be removed. Practically, we can notice that it comes 
from the exponential decay of the function b k , b k (0,L) ~ exp(— kL). The final expression for the 
path-integral after removing this artificial pressure is thus 



(2) 



exp 



27r6 fc (0,L) 



\ [« f fc (o,i) 



+ 



C + £cxtfc/47T 



-b k {o,Ly 



(A24) 



Appendix B: Fluctuations correlation function 

We show here how to compute the fe-mode correlation function Gk(z, z') and dd'Gk(z, z') needed 
in (|44| . 

First, we can derive that 

_ (ext,k\K i , MF (0,z;k)9K, l p MF (z,z';k)eK i , MF (z',L;k)\ext,k) 
h[Z,Z> (ext,fc|^ MF (0,L;fc)|ext,fc) ' 1 } 

which can be evaluated explicitly with (|A12I) and (|A8I) and leads to 

G h {,, z) = (aH0, z) - c+ee ££ > z) + [z -> L - z]j , (B2) 

where we took advantage of the symmetry z — > L ~ z of the system. 
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The term dd'Gk(z, z) can then be obtained as the limit 

dd'G k (z, z) = lim(2C)- 2 [G fc (z + £ z + C) - 2G k (z - (, z + () + G k (z - (, z - ()}. (B3) 

We can again compute this expression analytically, because we have an exact expression for 
K^, MF (z,z + (;k) when ( is small. As expected, it diverges. Again, we have to take the size 
R of the ions into account. To be consistent with the cut-off introduced for pressure, we will reg- 
ularize the divergence in Fourier space. The point is that locally, the correlator can be computed 
analytically, and that this local form contains the divergence; thus we will be able to cut it off in 
the analytic expression. 

In this purpose, we decompose the correlator in the following way 

G k = Gl um - G£ c + G£ c , (B4) 

where G£ um is the correlator computed using the numerical values for the Pauli-van Vleck functions 
(dealing with it in Fourier space would need a numerical Fourier transform) and Gj^ c is the " local" 
correlator, computed assuming that A(z) and u> k (z) are constant around z. Thus the term GjJ um — 
Gfc c will not lead to divergences when we will compute its second derivative, and we will be able 
to write the term G^ c in Fourier space and cut its divergence off easily. Explicitly, 

r loc, A CX P (-U fc (z)|z-Z , |) 

° k (Z ' Z)= 2A(z) Uh (z) ' (B5) 



its expression in Fourier space is 



G l g c (q) = -. . . „ r ^r, (B6) 

A(z)(q 2 +w fe (z) 2 )' 



so its contribution to the second derivative will be 



f 



nA(z) 



<?max(fc) - £*>fe(z) arctan 



(fe) 

LOk(z) 



(B8) 



where the superscript "cut" means that the cut-off has been applied, and the onc-dimensionnal 
cut-off depends on the mode since the three-dimensionnal wave- vector is constrained: 

?max(fc) 2 + k 2 = ££ ax . (B9) 

Now, we turn to the evaluation of dd'G k Alm (z, z). The correlators invoked in the above expression 
can be derived from the quantity 

Z{a,a',b) = J exp ^0' 2 + W0^ d9d9' = ^_ ^ , (BIO) 



where 



a = c + a!°(z-(,z + (), (Bff) 
a' = c' + a£(z-C,z + Q, (B12) 
b = b k (z-(,z + (), (B13) 



25 



with c = af(0,z - C) ~ enjggj^ and c' = of (0,L - * - C) - ^Sl^j^y 

We note that the terms appearing in a, a' and b can be computed reversing the composition rules 
(|A131IA15I) . giving 

affz-Cz + Q = -7- 6fc( °' Z ~ff 2 -a f fe (0,z-C), (B14) 



fr fc (0,£-z~C) 2 
a?(0,Z-2-C) -af(0,£-^ + C) 



a?(«-£z + = fc , n F _ ZkTr : r — --af(0,L-z-(), (B15) 



b k (z -C,z + C) = y^f^l K fc (°' 2 - ~ - C, * + 0] • (B16) 

(B17) 

where we used a k (z — (, z + () = a k (L — z — (,L — z + (). Now 

oi-e-c-o--^-^. (bis) 

0! -(, + C,, + 0--2^- S? !^. (B19) 

GT^-C.i + O = 2^2= » . (B20) 

ao aa' — tr 

We will thus write 

where [99'] £ means that the second derivative is computed with a small step £. To this quantity, 
we have to deduce 



[dd' ]c GtM= 1 -Z { ^if ) - (B22) 



Finally, we have 

99'G fe (z,z) = [dd'] ( Gl am (z,z) ~ [99']c G L° c ( 2 ^) + 99 ' G fe Ut (2^), (B23) 
that should not depend on the discretisation step £ as soon as it is small enough. 



